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ABSTRACT 

A perturber may excite a coherent mode in a star cluster or galaxy. If 
the stellar system is stable, it is commonly assumed that such a mode will 
be strongly damped and therefore of little practical consequence other than 
redistributing momentum and energy deposited by the perturber. This paper 
demonstrates that this assumption is false; weakly damped modes exist and 
may persist long enough to have observable consequences. 

To do this, a method for investigating the dispersion relation for spherical 
stellar systems and for locating weakly damped modes in particular is developed 
and applied to King models of varying concentration. This leads to the following 
remarkable result: King models exhibit very weakly damped m = 1 modes 
over a wide range of concentration (0.67 < c < 1.5 have been examined). The 
predicted damping time is tens to hundreds of crossing times. This mode causes 
the peak density to shift from and slowly revolve about the initial center. The 
existence of the mode is supported by n-body simulation. 

Higher order modes and possible astronomical consequences are discussed. 
Weakly damped modes, for example, may provide a natural explanation for 
observed discrepancies between density and kinematic centers in galaxies, the 
location of velocity cusps due to massive black holes, and m = 1 disturbances of 
disks embedded in massive halos. Gravitational shocking may excite the m = 1 
mode in globular clusters, which could modify their subsequent evolution and 
displace the positions of exotic remnants. 



Subject headings: galaxies: kinematics and dynamics, interactions, nuclei- 
globular clusters: general 
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1. Introduction 



A stellar system responds to a disturbance by both phase-mixing, an incoherent 
response, and by collective excitation of discrete modes, a coherent response. Since galaxy 
and star cluster models must be stable to be useful, growing mode^ have been studied 
extensively both analytically in the form of stability criteria and using n-body simulation 
(see Binney & Tremaine 1987, hereafter BT, for a review). The most well-known of these 
is the bar-forming mode which significantly reorganizes its host galaxy. In general, damped 
modes will also be excited by a disturbance and transport momentum globally but their 
amplitude decays with time. If the decay time is short, these modes will not be perceptible. 
However, if a mode is weakly damped, it may persist for many dynamical times, even 
though the system is stable. For example, a galaxy "fly-by" could excite such a mode which 
might persist long after the encounter even if the system is not strictly unstable. Several 
authors, in particular Miller and Smith (Miller 1992, Miller & Smith 1993), have observed 
oscillatory modes in simulations. Mathur (1990) formally demonstrated the existence of 
oscillating solutions in both one-dimensional and spherical systems. In this paper, we 
explore the existence and implication of weakly damped modes in spherical stellar systems. 

The existence of damped modes is not in question. Collective damping in a 
homogeneous plasma. Landau damping, is well-known. The simultaneous solution of 
the Vlasov equation and the Poisson equation yield a dielectric response function which 
describes the reaction of the plasma to an imposed disturbance. Even in the presence 
of no external disturbance, a stochastic space-charge field may give rise to spontaneous 
density fluctuations resulting in Landau modes. The zeros of the dielectric function are the 
conditions for collective solutions. This use of the response function is sometimes known as 
the dispersion relation (e.g. Ichimaru 1973) since the it determines the frequency-wavelength 
relation for electromagnetic waves in a medium.0 The relationship between the full modal 
spectrum and the initial value for the homogeneous plasma has been discussed by van 
Kampen (1955) and the same physical principles apply here. 

Damped modes are difficult to study in arbitrary stellar systems for two reasons. 
First, the dispersion relation is analytically intractable in nearly all but the homogeneous 



^For our purposes, a mode is a well-defined pattern with a possibly complex frequency, u. If Im{uj) > 0, 
the initial equilibrium is unstable and the mode initially grows exponentially. If Imiuj) < 0, the equilibrium 
is stable and the modes are oscillatory or damped. 

^The analogy does not carry over the stellar dynamical problem but the term dispersion relation is still 
used. 



-3- 



cases. Second, although one can study damped modes in simulations, one must be able to 
discriminate the mode from numerical artifacts. In this paper, I will develop a numerical 
procedure for evaluating the dispersion relation for spherical stellar systems and locating 
damped modes (§2). Armed with a prediction for the shape and frequency of a damped 
mode, we may initialize and efficiently search for it in n-body simulations. Spherical models 
are considered because of the theoretical simplicity, because they represent a wide variety 
of astrophysical scenarios and to relate to the large body of hterature on spherical stellar 
systems. In §3, the technique is applied to King models of varying concentration. King 
models are known to be stable, for example, by Antonov's criterion (Antonov 1962, see also 
BT). Nonetheless, we will see that these models all have very weakly damped I — m — 1 
modes, with damping times longer than ~ 20 half-mass crossing times. This is contrary to 
the widely expressed belief that a disturbance in a stable system will phase mix in several 
crossing times. As a check, wc realize the I = m = 1 mode and follow its evolution in a 
set of n-body simulations. The results of the simulation (§4) show that the expected mode 
persists with no obvious decay. The paper concludes with a discussion of astronomical 
applications (§5). 



2. Method 

For coUisionless systems, the dispersion relation describes the solution to an initial 
value problem as the result of a Laplace transform; each component has time- dependence 
of the form exp{—iujt), where the frequency oo may be complex. Although the dispersion 
relation for infinite homogeneous stellar systems has been derived (Ikeuchi et al. 1974, BT), 
most stellar systems are finite and inhomogeneous. In a stellar sphere in particular, orbits 
are quasiperiodic and the repetitive effect of a disturbance at the same point in an orbit's 
phase can have a large effect (resonance) . Periodicity may be introduced in a homogeneous 
medium by carving out a cube and connecting opposite faces; topologically, this is a 3-torus. 
The resulting so-called periodic cube has been used by several authors (Barnes et al. 1986, 
Weinberg 1993) to take the quasiperiodicity into account without sacrificing the theoretical 
simplicity of spatial homogeneity and will be used below as an example. 

Inhomogeneous models of galaxies and clusters are considerably more difficult to treat 
analytically because the coUisionless Boltzmann and Poisson equations are separable in 
different bases. By Jean's theorem, the solution to the coUisionless Boltzmann equation 
may be written conveniently in terms of its integrals of motion. In the spherical case, these 
are energy and angular momentum or radial and angular actions. Since the unperturbed 
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Hamiltonian is then independent of associated angle-like coordinate variables, a linear 
perturbation to this equilibrium may be written as a Fourier series in the coordinate 
angles. The Poisson equation is naturally expanded in a different basis: spherical hamonics. 
However, by expanding the self-consistent gravitational potential in an harmonic series 
which satisfies the Poisson equation by construction, the simultaneous linearized solution to 
the whole system may be reduced to a matrix equation. In particular, let a^, bj be two sets of 
expansion coefficients. Then the response bj to a perturbation aj is given by bi = Mij{u)aj 
where Mij is the response matrix and describes the dynamical response of the stellar sphere 
to that part of perturbation with time-dependence exp(— icut). A self-consistent response 
has Oj = bj (cf. Appendix B). This formalism, sometimes called the matrix method, has 
had varied applications beginning with Kalnajs (1977) who was interested in the modes of 
stellar disks. Polyachenko and Shukhman (1981) first adapted the method to the study 
a spherical system (see also Fridman and Polyachenko 1984, Appendix) and the method 
was later employed by both Palmer and Papaloizou (1987) in the study of the radial orbit 
instability and by Bertin and Pegoraro (1989) to study the instability of a family of models 
proposed by Bertin and Stiavelli (1984). Weinberg (1989) used the matrix formulation 
to study the response of a spherical galaxy to an encounter with a dwarf companion and 
Saha (1991) and Weinberg (1991, hereafter Paper I) investigated the stability of anisotropic 
galaxian models. Using this approach, the dispersion relation is the condition that the 
self-consistent response have a non-trivial solution: det[5ij — Mij^u)] = 0. The general 
dispersion relation can be calculated for a spherical model with an arbitrary phase-space 
distribution function, although only isotropic distributions will be discussed here. 



2.1. Analytic continuation of the dispersion relation 



Just as in the case of the plasma dispersion relation, the derivation in Paper I 
explicitly shows that the dispersion relation as written (see Appendix B) is only valid in the 
upper-half u plane. Causality in the initial value problem demands that various integrands 
be analytically continued to the lower-half plane. This is readily done in the plasma case 
(cf. Ichimaru 1973) where the natural coordinates for both the Boltzmann and Poisson 
equations are cartesian. This is problematic in the case of the sphere, whose dispersion 
relation is given by equations (p5|), ( |B6|) , and (pTj). Since equation (p5| ) contains the term 
l/[a; — n ■ 0(I)], explicit analytic continuation requires analytic continuation of orbital 
frequencies and actions. To sidestep this daunting task, I propose to evaluate the dispersion 
relation in the upper-half plane and perform the analytic continuation numerically by 
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approximating it by a function which may be easily continued. Given that poles are 
expected because of the form of equation ([B5|), rational functions are a natural choice. 

Before applying this to cases of interest, we would like some assurance that this 
procedure for analytic continuation will work reliably. To do this, I evaluated the dispersion 
relation for the periodic cube (eq. |Al| ), both analytically and numerically. The simple form 
of the velocity integrals in equation allow explicit analytic continuation. Analytically 
generated contours for the fundamental harmonic in the dispersion relation, 11)10,0(1^)1; 
are shown in Figure |I|. The zeros are at the centers of the concentric contours and the 
location of the zeros indicate the type of growth; if the zeros are in the lower (upper) 
half-plane, disturbances are damped (growing). The model depicted is stable. The most 
weakly damped mode has Re{uj) = and and therefore damps in situ. The other modes 
both damp and propagate with group velocity Re{uj). As the stellar velocity dispersion is 
decreased, the solution with Re{uj) = moves toward the real axis from below and will 
become the Jeans instability for sufficiently small velocity dispersion. Modes in spherical 
systems, discussed below, have spatial distributions proportional to spherical harmonics 
and Re{uj) indicates a pattern speed rather than a group velocity. 

I then evaluated Pi,o,o(^^) at 20 points for —40 < Re{uj) < 40 for each of 
Im{(jj) = 0,0.5, 1.0 and found the corresponding rational function fit (cf. Appendix p] ). 
Contours of the approximating rational function are shown in Figure ^. Comparing Figures 
m and H we see that the least damped modes are determined accurately. I checked the 
robustness of the method to the choice of evaluations by varying the number and placement 
of points included in the rational function fit. The three most weakly damped modes are 
accurately determined even if any two of the "rows" of constant Im{(jj) are eliminated. The 
most weakly damped mode remains accurately determined even for 10 points along the real 
cij-axis from —40 < Re{uj) < 40. The behavior of T>i q^q near the real axis is dominated 
by the true weakly damped modes, allowing zeros to appear at Im{uj) ^ —28 in Figure 
^. This is an expected and generic feature; zeros on the lower-half plane will only be 
well-approximated by the rational function if they affect the V on the upper-half plane. In 
general, we can not expect accurate determination of all features on the lower half-plane 
from a finite grid and, in particular, the upper half-plane will place little constraint on the 
rational function extrapolation on the lower half-plane if more weakly damped modes at 
similar values of Re{u}) exist. 



2.2. Application to the dispersion relation for stellar spheres 



- 6 - 



In the next section, we investigate the low-order damped modes for isotropic King 
models with dimensionless central potential Wq = {Et~ Eo) / a"^ = 3, 5, 6, 7 and concentration 
c = log^Q(r(/rc) = 0.67, 1.0, 1.3, 1.5 respectively where and are the core and tidal 
radii (see King 1966 and BT for more details). Linear perturbation theory allows each 
harmonic order, oc Yim, to be considered separately. However, because the sphere has no 
symmetry axis, all terms m are degenerate for a given /. In realizing specific modes below, 
we will choose I = m. Units are chosen so that M = G = 1 and the total gravitational 
potential W = —1/4. In most cases, we use a grid in the upper-half cj-plane with 
Im{uj) = 0.01, 0.05, 0.1 and 10 to 40 evenly spaced points in —1 < Re{uj) < 1. In some cases 
the domain is extended to —5 < Re{uj) < 5 to check for "fast" pattern speeds but none were 
found in any cases explored here. The evaluations of T>{uj) (cf. eq. [B7| ) are fit to a rational 
function and zeros in the lower-half plane are noted. The robustness of the zeros to the 
input grid is checked by successively truncating the grid in | Re{uj)\, halving the resolution 
of the grid and leaving out "rows" in Im{uj). In all cases, the existence of the zeros and 
general topology of the function did not change although the position of particular zeros 
varied by as much as 10% in a few cases. 

The elements of the response matrix Mij{uj) (cf. eq. p5|) require a two-dimensional 
integral and a two-dimensional sum in li and l2- The sum is only two dimensional since 
Mij{(jj) is independent of /s = m given /. The sum in li must be truncated at a chosen 
value, li-max and the sum in I2 has I/2I < ^- Twenty point grids for the integral evaluation 
using rational functions were found to be adequate; large numbers of grid points 40) 
cause significant truncation error (see App. D). The rational function technique gives 
surprisingly good results with a small number of integrand evaluations. Standard lower 
order integration techniques (e.g. Simpson's rule) ultimately give better results but only 
after an order of magnitude more evaluations. The value of h^max was successively increased 
to obtain convergence. The appropriate value of h^max depends on I; h^max was chosen to 
be 6 and 10 for / = 1 and 2, respectively. Finally the order of the n x n matrix Mij was 
varied to obtain convergence in V{uj). The required order increases with Re{uj); n = 30 
gave adequate convergence for | Re{uj)\ ^ 1. The resulting error in Vi^uj) is estimated to be 
0(10-2). 



3. Properties of damped modes 

This section describes the application of the methods from §2 in determining the 
damped modes of King models with Wq = 3, 5, 6, 7. The complex frequencies of the most 
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Table 1: Model parameters and modal frequencies 



Model 




Order 


Frequency 


Shift 


Example 


King Wq 


c 


/ 


Re{uj) 


Im{ijj) 


r shift/ 


Cluster name 


3 


0.67 


1 


0.010 


-0.0016 


0.41 


NGC 6496 






2 


0.10 


-0.076 






5 


1.0 


1 


0.034 


-0.0011 


0.40 


NGC 288 






2 


0.0 


-0.11 






6 


1.26 


1 


0.042 


-0.0065 


0.42 


oj Cen 


7 


1.5 


1 


0.064 


-0.025 


0.54 


M13 






2 


0.0 


-0.20 







weakly damped modes are shown in Table 1. To summarize, one finds the / = 1 modes are 
very weakly damped for all Wq studied with the damping proportional to concentration. 
The existence of such modes imply that low-order excitations of stable galaxies and clusters 
may persist for many crossing times and scaled to physical units, perhaps as long as the age 
of the galaxy itself. Very weakly damped / = 2 modes were found for the Wq = 3,5 cases 
and the same mode is relatively more strongly damped for Wq = 7. No weakly-damped 
/ = 3 modes were found. A detailed discussion of these trends and an examination of the 
modal shape is given below. 



3.1. Frequencies of lowest order modes 



For comparison. Figure |^ plots the zeros of the dispersion relation for / = 1 damped 
modes in Wq = 3, 5, 6 and 7 King models. For example, the I = 1 harmonic for the Wq = 5 
King model has a very weakly-damped mode at = (±0.034, —0.0011). The tu-plane was 
searched for —1 < Re{uj) , Im{uj) < 1 and no other modes were indicated. Changes in the 
input grid of T>{ijj) change the location of Im{uj) for the modes by (9(10"'^). The errorbar 
in upper left corner of Figure ^ gives the estimated precision in u. Both the Wq = 3 and 5 
have nearly zero damping within the resolution of the technique; the value Im{uj) = can 
not be ruled out for these cases. The analogous modes for the Wq = 6, 7 models are clearly 
damped but still only weakly. 

The I = 2 modes exhibit a similar trend: the rate of damping increases with 
concentration. For the Wq = 5, 7 cases, the pattern speed is zero but it is nonzero for the 
lowest concentration model. 
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Grids for the I — 3 harmonic were constructed for models Wq — 3,5,7 but no weakly 
damped modes with | Re{u!) \ ^ 0.3 and Im{uj) ^ —0.3 were found, suggesting that the I — 3 
disturbances are moderately to strongly damped. 
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3.2. Shape of lowest order modes 

The / = 1 weakly-damped mode for the Wq = 5 King model discussed in is shown 
in Figure The procedure for finding the eigenfunction is discussed in Appendix ^ The 
peak of the disturbance occurs at ~ 75% of the core radius of the background model. The 
half-mass radius is ri/2 ~ 1-6. At larger radii, the density contours remain concentric. 
However, the shape of an m = 1 mode depends on the expansion center; therefore one 
would observe this mode as a coherent m = 1 shift at large galactocentric radius by fixing 
the expansion center on the density center. 

The / = 1 mode manifests itself as a shift of the central density. The shift is given 
by the linear relation r shift/ fc ~ 0.4e, where e is the ratio of perturbed density at peak 
relative to background and e ^ 0.3 to have positive density everywhere. As an example, the 
solid contours in Figure ^ show the combined background and perturbation with the near 
maximum amplitude of 30% at peak. The dotted contours show the unperturbed model 
at the same levels for comparison. The shifted center revolves about the original center 
at angular frequency Re{uj). The center-of-mass shift for other concentrations is given in 
Table 1. The value r shift /^c^ is similar for all concentrations. 

The magnitude of the perturbed phase-space distribution averaged over angles is shown 
in Figure ^ The energies of perturbed orbits are small but do not include the most bound 
particles and the peak toward large at constant E indicates that the perturbed orbits 
tend to be more circular on average. 

The I = m = 2 mode (Fig. ^ appears as a bar which peaks at a fraction of the core 
radius, at roughly the same radial scale as the I = 1 mode. Figure || shows the mode at 30% 
amplitude combined with the Wq = 5 King model background density. Since Re{uj) = 0, 
the pattern will have a fixed orientation and the bar distortion will decay with time. The 
phase space perturbation is similar in extent to that shown in Figure ^ but with a larger 
contribution of eccentric orbits. 



4. Simulation 

Although the dispersion relation indicates the existence of a weakly damped / = 1 
mode, can such a mode be excited at significant amplitude? 

In partial answer, we may attempt to realize the mode in an n-body simulation. Success 
would also serve as an independent check of the linearized solution. Unfortunately, such a 
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simulation is fraught with serious systematic problems. First, most simulation schemes do 
not explicitly conserve momentum. Since an / = 1 perturbation shifts the density center, 
numerical artifacts may not be distinguishable from the mode itself. To minimize possible 
artifacts, we use a direct-summation force calculation which does conserve momentum. 
Second, minimization of two-body relaxation requires softening at, roughly, the mean 
interparticle spacing. However, the feature we wish to resolve is only a factor of a few larger 
than this scale for ~ 10000 bodies. Therefore the softening itself is likely to change the 
dispersion relation and this will require that the dependence on softening be investigated. 
Third, even with softening, relaxation significantly changes the background on a shorter 
timescale than the predicted damping, suggesting that relaxation in the simulation may 
affect the mode. Finally, although these difficulties could be addressed with larger n, the 

FLOP scaling of the direct summation scheme practically limited the simulations to 
n = 10000.0 

Because the mode is computed in action-angle variables, the perturbed n-body phase 
space is easiest to realize in these variables. The details are described in Appendix |E] The 
mode described in §3.2 is realized with a peak level of 20% of background. The softening is 
chosen to be the mean interparticle spacing at the half-mass radius. Because of both the 
imposed perturbation and modification of the interparticle force by softening, the initial 
conditions are not in equilibrium. The system appears to settle at t = 25 as indicated by 
a constant value of —2T/W ~ 1 (cf. Fig. For comparsion, the radial crossing time at 
the half-mass (~ 1.6) is roughly 4 time units. The —2T/W profile looks similar for initial 
conditions without an initial perturbation, suggesting that softening causes the poor initial 
equihbrium state. 

The evolution of the introduced mode may be diagnosed using the same harmonic 
decomposition used to evaluate the dispersion relation (Appendix B). The original expansion 
center is preserved in the center-of-mass frame since the direct force scheme explicitly 
conserves momentum. The azimuthal phase of the I = m = 1 distortion for several of 
the higher-order radial wavefunctions are shown in Figure 0. The low-order radial terms 
are less sensitive to the mode which has a scale length of ~ 1/10 the system radius. The 
real part of the eigenfrequency is the pattern speed for an m = 1 distortion. The dashed 
line indicates the slope predicted by Re{uj). The agreement is quite good. However, the 
evolution due to relaxation removes sensitivity to changes in amplitude of the mode which 
is predicted to have a longer timescale. 



•^Computations were performed on a Sun Microsystems Sparc 10/41 workstation which required 3.4 x 10^ 
CPU seconds per step. 
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The amplitude of a disturbance may be computed from the total gravitational potential 
energy in the disturbance. From Appendix B, it follows that this energy for a particular 
harmonic is Wi^ — —2t:GY^^=\' where a'*" are the expansion coefficients. For the 

current simulation, the energy per order /, Wi = J2m ^im, converges with increasing radial 
wavenumber Umax at fixed time, confirming that the disturbance is global. In addition, Wi 
is roughly constant over the course the simulation. More precisely, Wi does appear to grow 
slightly with time but this can not be distinguished from the effects of relaxation which 
causes the central potential to deepen. An additional simulation with a 7% rather than 
20% amplitude perturbation exhibited similar behavior, both in phase and amplitude. 

To understand the effects of softening, simulations with softening length of both 
half and double the half-mass interparticle spacing were performed with the same initial 
conditions. The same general trends were seen: the introduced I — 1 mode persisted with 
the expected pattern speed. 

However, in all three runs, there were "bursts" of the I — 2 mode, which might be 
expected given its long damping time (cf. Table 1). The I — 2 power was larger for smaller 
softening, possibly suggesting that either or both relaxation and softening are changing the 
modal structure. In addition, for t ^ 30, relaxation causes the central potential to deepen 
and at the same time, the relative contribution of the / = 2 disturbance grows and appears 
to saturate without damping, roughly at the level of the initial / = 1 perturbation. For the 
case with 7% perturbation intially, the I = 1 mode continues at the predicted pattern for 
t ^ 30 even though the / = 2 distortion is evident. For the 20% case, the growth of the 
/ = 2 distortion appears to significantly change the / = 1 mode. The importance of these 
effects (especially the possible implication for bar growth) remains to be investigated. 

Finally, a set of simulations with no perturbation a priori were used to check the 
significance of the observed disturbance by comparing the gravitational potential energy in 
/ = 1 mode. The stochastically excited / = 1 components contained ^ 3% of the energy 
in the introduced mode described above. Nevertheless, it is worth noting that a coherent 
pattern was seen in these control runs with the expected pattern speed even though the 
peak density contrast was only 1% to 3%. 



5. Discussion and Summary 

The numerical evaluation of the dispersion relation for King models predicts the 
existence of low-order weakly damped modes. In particular, the 1 = 1 (m = 1) mode is very 
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weakly damped with small pattern speed over a range of central concentrations. Similar 
results might be expected for other spherical systems. 

This mode causes the density center to shift and rotate about the original center. Using 
the linear theory, the shift is given by the linear relation r/rc ~ 0.4e, where e is the ratio of 
perturbed density at peak relative to the background and must be ^ 30% of background at 
peak. This mode has been realized and observed in a simulation of a Wq = 5 King model. A 
strong excitation of this mode results in a shift of the density center of roughly 10% of the 
core radius as predicted and such levels were maintained for the duration of the simulation. 
However, intrinsic limitations and features such as softening and relaxation modify the 
evolution of and limit temporal sensitivity to weakly damped modes. This, in part, explains 
both the reason they were not discovered by n-body simulation and the general belief that 
excitations of stable equilibria vanish in several crossing times. 

Because the low-order damped modes are simple, they should result from general 
perturbations. For example, the differential acceleration due to a passing galaxy may cause 
a significant m — 1 perturbation which might mix quickly except for the weakly damped 
modes. A companion on an eccentric orbit may have similar effects. The induced distortion 
due to a companion on circular orbit will be dominated by the quadrupole {I — 2) moment; 
the dipole force is canceled by the barycentric motion (Weinberg 1989). 

In the halo of a halo-dominated galaxy, such a mode could result in the appearance of 
an m = 1 distortion or position-shift of the luminous disk. Baldwin, Lynden-Bell and Sancisi 
(1980) have pointed out that a number of nearby galaxies have lopsided HI distributions 

relative to their optical centers and that a large fraction of spiral galaxies show this effect. 
Damped modes may be partly responsible for the ubiquity of such offsets, especially in the 
more moderate cases. Using the Wq = 5 model as an example, the excitation of a 10% 
perturbation in a halo with core radius of 15 kpc by, say, a companion or "fiy-by" would 
result a central density offset of ~ 600 pc which would persist long after the encounter. In 
the bulge or spheroid of a disk galaxy, the mode would shift the density center from the 
kinematic center as determined, say, by the large-scale HI rotation curve. Using the same 
example, for a bulge-spheroid with a 1 kpc core radius with e = 0.01, 0.1 the offset would 
be ~ 4, 40 pc. In addition, since the mode is slowly varying, one might expect a massive 
central black hole, M ^ 10^ M© to remain within the potential well of the shifted density 
center since the equipartition radius would be rcm^/M},h for mean stellar mass, m*. 

Similarly, a time-dependent external force or gravitational shock differentially 
accelerates globular clusters, possibily exciting weakly-damped modes. Clusters on eccentric 
orbits will suffer gravitational shocking due to the halo and bulge and all clusters with 
mean galactocentric radii within 20 kpc will suffer shocking due to the disk. In addition. 
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there are relaxation-driven internal momentum sources such as binary recoil which could in 
principle produce an m = 1 disturbance (Makino & Siigimoto 1987). Either way, an offset 
core may change the subsequent dynamical evolution by decreasing the thermal contact 
of the core and the population of eccentric halo orbits. This may also have implications 
for gravothcrmal oscillation. In addition, we may have to reexamine the intepretation of 
off-centered compact objects if any of the m = 1 excitations keep the density central offset, 
on average, from the geometric center. On the other hand, the physical picture is further 
complicated by two-body relaxation whicy may lead to more rapid damping of these modes. 
Also, the largest relative central density shifts and damping times occur for the lowest 
concentrations where we expect rates of relaxation and subsequent production of exotic 
objects to be low. Additional work will be necessary to understand these competing trends 
and the extent to which these modes may be excited in globular clusters. 

I thank Susan Kleinmann and David van Blerkom for comments on the manuscript 
and Dave Chernoff and Ira Wasserman for discussion. This work was supported in part by 
NASA grant NAGW-2224. 
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A. Dispersion relation for a stellar cube 



As discussed in §2, the dispersion relation is the condition for simultaneous solution of 
the collisionless Boltzmann and Poisson equations in response to a perturbation. In the case 
of a homogeneous cube with periodic boundary conditions, the two equations are separable 
in the same coordinate system, cartesian, which leads to an analytic solution. Following 
Barnes et al. (1986), the dispersion relation may be written: 

= 1 + — ^An, (Al) 

where n is a three vector of integers. 



3 k-dfo/dv 



An = t d'v^^ , (A2) 

M is the total mass of stars in the cube with side length L, k is the wave vector defined by 
k = 27rn/L, and fo is the background phase-space distribution function. The eigenfunctions 
are proportional to exp(ik ■ x). Barnes et al. (1986) derive equation ([Al|) by explicitly 
introducing a vanishing small growth term. For an alternative derivation using the Laplace 
transform see Weinberg (1993). 



B. Dispersion relation for the stellar sphere 



Any disturbance in the sphere may be represented as a spherical harmonic expansion 
with an appropriate set of orthogonal radial wavef unctions. We choose a biorthogonal 
potential-density pair as described in Paper I. The pair (m-"^, d-"^) is constructed to satisfy 
Poisson's equation, V^m-"^ = inGd''^, and to form a complete set of functions with the 
scalar the product 

— Urr\^*d["' = {^ '^' = { (Bl) 



47rG J « J [ g otherwise. 

The perturbed density and potential then have the following expansions: 



$1 = E>^^-(^></')E«fW<W' (B2) 

Im j 



Im 



Paper I shows that the coupled system of collisionless Boltzmann and Poisson equations 
may be solved by reducing their simultaneous solution to a matrix equation using such 
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an expansion above. The solution may then be written as the non-trivial solution to the 
following linear equation 



-a'r = M,^{u)~af, (B4) 

where summations over like indices are implied, the '~' indicates a Laplace-transformed 
quantity, 

{2-nf f f dEdJJ ^ 2 dfo 1 



-n 



A-nG J J VLi{E,J)^ 21 + 1 dluj-u-Q. 



F,Jvr/2,0)ri^/j::,;(I)H^^l:f. (I), (B5) 



and 



Wlli^{h,h) = - r dw^cos[hw^~h{'iP-W2)]uf{T). (B6) 
7r JO 



In the above equations, Ij are the actions, Wj are the angles, and n = {h,l2, h) is a vector 
of integers. The actions were chosen by solution of the Hamiltonian-Jacobi equation; Ji is 
the radial action, I2 is the total angular momentum and I3 is the z-projection of the angular 
momentum. The quantity — W2, the difference between the position angle of a star in 
its orbital plane and its mean azimuthal angle, depends only on wi, the angle describing 
the radial phase. The frequencies associated with the conjugate angles w are defined by 
ft = dH/dl. Qi is the radial frequency and reduces to the usual epicyclic frequency for 
nearly circular orbits. Q2 is the mean azimuthal frequency. ^3 = and the corresponding 
angle is the azimuth of the ascending node (see Paper I or Tremaine and Weinberg 
1984, for details). It follows that M*j{lj) = Mij{—u*) and this property is exploited in the 
computations described in §3 to reduce the number of numerical evaluations. Since most 
researchers quote phase-space distribution functions in energy and total angular momentum 
f{E, J), I have transformed the integration in the definition of M(ci;) in equation ( P5| ) to E 
and J variables. The factor l/[u; — n ■ 0,(1)] in the integrand of equation ( |B4| ) may cause 



discontinuities in the E and J integration. This is efficiently addressed using the rational 
function techniques discussed in Appendix |D] 

Equation ( |B4|) only has a nontrivial solution if 

V{iu) = det{l - M{lu)} = 0. (B7) 

As in plasma theory, we refer to the function T>{uj) as a dispersion relation. In general, 
uu will be complex. Recall that the coefficients dj which appear in equation ( [B4| ) are 
Laplace transformed and describe the response of the system to a perturbation with a 
time-dependence of the form exp{—iujt). Because equation ([B4[) follows from a Laplace 



transform, it is only valid in its present form for 00 in the upper-half complex plane. The 
dispersion relation must be analytically continued to Im{uj) < to find damped modes. 
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The response coefficients a''" which solve equation ( |B4|) for a particular eigenfrequency tu, 
describes a mode of the system. 

Because a sphere has no unique symmetry axis, a mode can not depend on an arbitrary 
choice of coordinate axes. Since a rotation causes m components to mix according to the 
rotation matrices, the mode itself must depend on I alone. In particular, note that since 
the functions (m'™", d'™) may be chosen independent of m (e.g. the spherical Bessel function 
used by Fridman and Polyachenko 1984, and Weinberg 1989), equation ( [B5| ) and thus the 
dispersion relation is independent of m. 



C. Eigenfunctions 

The eigenfunction corresponding to the frequency at which T>{u)) = may be obtained 
by simultaneous solving the n — 1 independent equations implied by the n x n matrix given 
by equation ( |B4| ) for the expansion coefficients aj. The density and potential perturbations 
follow directly using equations ( |B2D and (PB|). Recall that the eigenvalue is found by 
rational function fit to a grid of T>{uj). The matrix elements at the eigenvalue are determined 
by rational function fit to the grid used to find the eigenvalue originally. The dispersion 
relation is not solved exactly by the new interpolated matrix Mij[uj) and is iteratively 
refined to find a new solution to Viuj) = 0. Reassuringly the new eigenvalue is found within 
5% in \u\ of the predicted value in all cases. 



D. Rational function techniques 

Rational functions are used for two tasks in the this paper. First, they are suitable 
estimators for the functional form of the dispersion relation. Using a discrete number of 
evaluations of the dispersion relation in the upper-half plane, the rational function may be 
trivially analytically continued to the lower-half plane and, in particular, zero locations may 
be straightforwardly determined. Second, rational functions are ideally suited to evaluating 
complex integrands with poles or near singularities. Since the integrands of equation ( P^ ) 
have "vanishing" denominators — denominators of the form 1/[uj — n ■ fi(I)] — one expects 
that rational function to approximate these integrands rather well. 

The general procedure is as follows. Suppose we are given n evaluations of a 
complex function f{zi). The standard techniques for rational function evaluation may be 
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straightforwardly extended to the complex domain. In particular, I adopted reciprocal 
differences along the main diagonal (Stoer & Bulirsch 1980) from which the coefficients of 
the numerator and denominator polynomials may be determined directly by recursion. Let 
the rational function be written as 

m = ^. (Di) 

where N and D are polynomials and let the order of be n and the order of D be d. The 
reciprocal difference routine gives either n = d or n = d+ l. Using synthetic division, 

R{z) = Q{z) + ^ ^ Q{z) + R{z) (D2) 

where the order of Q and will be and n — 1 if n = d or 1 and n — 2iin = d+ l. Now 
since i?(oo) = 0, one can show that 

k 

R{z) = J2S.{z) (D3) 

i=l 

where the sum is over the k zeros in D{z), Zi, and Si{z) is the principal part of R{z) at Zi. 
Following Henrici (1974), we may write the principal as follows: 

rrii 

S,iz)=y2a,,jiz~Zir^ (D4) 
i=i 

where rrii is the multiplicity of the zero at Zi. 

An efficient method for computing the principle parts has been outlined by Henrici 
(1974); I will sketch his arguments here for completeness. Since the coefficients of R and D 
are available, the poles of D may be derived numerically using deflation and R{z) may be 
written in the form 

^^\z-ZiJ di[z) 

By construction n{z)/di{z) is analytic at Zi and may be expanded in non-negative powers 
of h = z — Zi and therefore 

oo 

Riz, + h) = ^Q,„/i"-'"\ (D6) 

n=0 

Comparing equation (D6) with equations (|D3| ) and (|D4|) , it follows from the uniqueness of 
Laurent expansion that j for j = 1, . . . , mj. The desired values of may be 

then straightforwardly obtained from the explicit expression of n{z) and di{z). 
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This approach is remarkably useful since any rational function may be integrated 
exactly, limited only by truncation errors and the accuracy of the initial polynomial 
coefficients. Therefore, once an expression in the form of equation (|D2|) is obtained, the 



rational function R maybe be trivially integrated along any segment of the real line. In 
practice, I have found that rrii is rarely larger than 1 and the computational labor is usually 
much smaller than the statement of general method might suggest. 

One unfortunate feature of the rational function scheme is that does not converge with 
increasing grid density due to truncation error. For example, imagine trying to approximate 
a lower order rational function with a very dense grid of points. The exact representation 
will have identical, and therefore canceling, roots in both the numerator and denominator. 
Since these will not cancel precisely in the numerical representation, a very dense grid 
may result in a noisy approximation. For zero location, incomplete cancellation is easily 
recognized. For integration, one must check the asymptotic behavior. 



E. N-body realization in action-angle variables 



The most common procedure for realizing an arbitrary but spherical phase-space 
distribution is the rejection method for psuedorandom variables (e.g. Press et al., 1992). 
The procedure is as follows. Let TZ he a. variate from the unit interval. If the phase-space 
volume is finite for the system, the radius may be computed from the cumulative mass 
distribution: M(r) = TZM{rmax)- The gravitationally potential (f){r) then fixes the 
maximum velocity. A trial velocity may then be chosen uniformly from the sphere of radius 
Vmax = \J'^[4>ijmax) — (t'ij')] f^om which the energy and angular momentum, E and J, may 
be chosen. Let fo{E, J) be the phase-space distribution function. The trial point is accepted 
if /o(-E', J)/ sup |/o(-E, J) I > IZ. Otherwise the entire procedure is repeated. 

The perturbed phase-space distribution, / = /o + e/i, is best realized in action-angle 
variables but otherwise the procedure is analogous. From Paper I (see also Appendix B), 
the perturbed mode may be written 

A(I, w) = X: n ■ ^ ^_^^.^ e-("— ^) E o^r^MJ^- (El) 

One chooses the three trial actions, or equivalently E, J, J^, from an isotropic distribution, 
fcomp{E), chosen to lie everywhere above /. The angle variables are chosen uniformly from 
the interval [0,27r]. The orbit corresponding to the trial phase-space point is then specified 
in the background model. The phase point is accepted if f(l,w)/fcomp{E) > IZ. This 
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procedure is efficient for a careful choice of fcomp- The quantity e must be chosen to be 
sufficiently small that / > everywhere. 
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Fig. 1. — Contours in magnitude of the analytically-generated dispersion relation for the 
stellar cube. 

Fig. 2. — Same as in Fig. |l| but for the numerically-generated dispersion relation for the 
stellar cube. 

Fig. 3. — Shows the location of the damped / = 1 modes for Wq = 3, 5, 6, 7 King models 
(labeled by Wq). The estimated precision is indicated by the errorbar at the upper left. 

Fig. 4. — Contours of the density perturbation for the I = m = 1 mode for the Wq = 5 
model on the x-y plane. Only the inner region of the model is shown; the tidal radius is 
Tt = 8.72. Overdensity (underdensity) is shown as a solid (dashed) line. 

Fig. 5. — Effect of the density perturbation shown in Fig. ^on the background model. The 
peak perturbation is 30% of the background. Note the shift of the central density peak. The 
core radius is at 0.81 in these units. The dotted contours show the unperturbed background 
at the same levels. 

Fig. 6. — Absolute value of the phase-space perturbation for the I = m = 1 mode for the 
Wq = 5 model as a function of energy E and = {J / JmaxY- 

Fig. 7. — Contours of the density perturbation for the I = m = 2 mode in the Wq = 5 
model on the x-y plane. 

Fig. 8. — Effect of the density perturbation shown in Fig. |^on the background model. The 
peak perturbation is 30% of the background. 

Fig. 9. — The virial quantity —2T/W (top panel) and the run of gravitational potential 
energy (bottom panel) shown clS 8b function of time. 

Fig. 10. — Position angle as a function of time for the I = m = 1 component in the n- 
body simulation described in §4. Three radial wavefunctions (orders n = 8, 9, 10) are shown 
separately (solid lines) along with the theoretical prediction from §3 (dashed line). 



